home *** CD-ROM | disk | FTP | other *** search
- Frederick W. Hegeman
- 1022 South St.
- Rapid City, SD 57701
- (605) 343-7014
-
- Code to accompany article on factorial-base math
- ARITHMETIC IN FACTORIAL-BASE
-
- /* Listing 1. */
- /* facbase.c -- arithmetic in factorial-base */
-
- #include <stdio.h>
-
- #ifndef TRUE
- #define TRUE 1
- #endif
- #ifndef FALSE
- #define FALSE 0
- #endif
-
- #define DEBUG TRUE /* include RPN calculator */
- #define PLUS FALSE /* sign values */
- #define MINUS TRUE
- #define MAXFBINTEGER 100 /* <= 180 for arrays of ints */
- #define MAXFBFRACTION 100 /* <= 180 for arrays of ints */
- #define HIGHGUARD MAXFBINTEGER + 1
- #define LOWGUARD MAXFBFRACTION + 1
- #define ARRAYSIZE 1 + HIGHGUARD + LOWGUARD
-
- /* two constants in the proper format */
- int zero[ARRAYSIZE] = { PLUS, 0 }; /* signed 0 */
- int one[ARRAYSIZE] = { PLUS, 1, 0 };
-
- /* a flag set by divide() while "estimating" */
- int nowarning = FALSE;
- /* assign the value of one array to a second array */
- assignto(dest, source)
- int *dest;
- int *source;
- {
- int i;
-
- for(i = 0; i < ARRAYSIZE; i++)
- dest[i] = source[i];
- }
-
- /* z = abs(x) */
- absolute(x, z)
- int *x;
- int *z;
- {
- if(z != x) /* if not changing in place */
- assignto(z, x); /* copy */
- z[0] = PLUS; /* force positive */
- }
-
- /* Assign unary negative of x to z */
- negative(x, z)
- int *x;
- int *z;
- {
- int sign;
-
- sign = (x[0] == PLUS) ? MINUS : PLUS;
- if(z != x) /*if not changing in place */
- assignto(z, x); /* copy */
- z[0] = sign; /* change sign */
- }
-
- /* compare 2 factorial-base numbers for |x| > |y| */
- int greaterthan(x, y)
- int *x;
- int *y;
- {
- int i, j;
-
- /*
- * check the integer part first,
- * largest subscript to smallest
- */
- for(i = MAXFBINTEGER; (x[i] == y[i]) && (i > 1); --i)
- ;
- if(i != 0)
- {
- if(x[i] > y[i])
- return(TRUE);
- if(x[i] < y[i])
- return(FALSE);
- }
- /*
- * The integer portions are equal,
- * continue with the fractional part,
- * smallest subscript to largest.
- */
- for(i = HIGHGUARD + 1, j = 1;
- (x[i] == y[i]) && (j < MAXFBFRACTION - 1);
- i++, j++)
- ;
- return((x[i] > y[i]) ? TRUE : FALSE);
- }
-
- /* Compare 2 factorial-base numbers for |x| < |y| */
- int lessthan(x, y)
- int *x;
- int *y;
- {
- int i, j;
-
- /*
- * check the integer part first,
- * largest subscript to smallest
- */
- for(i = MAXFBINTEGER; (x[i] == y[i]) && (i > 1); --i)
- ;
- if(i != 0)
- {
- if(x[i] < y[i])
- return(TRUE);
- if(x[i] > y[i])
- return(FALSE);
- }
- /*
- * The integer portions are equal,
- * continue with the fractional part,
- * smallest subscript to largest.
- */
- for(i = HIGHGUARD + 1, j = 1;
- (x[i] == y[i]) && (j < MAXFBFRACTION - 1);
- i++, j++)
- ;
- return((x[i] < y[i]) ? TRUE : FALSE);
- }
-
- /* Compare 2 factorial-base numbers for |x| == |y| */
- int equalto(x, y)
- int *x;
- int *y;
- {
- int i;
-
- /*
- * Check the whole array except the sign
- * -- order doesn't matter.
- */
- for(i = ARRAYSIZE-1; (x[i] == y[i]) && (i > 1); --i)
- ;
- return((x[i] == y[i]) ? TRUE : FALSE);
- }
-
- /*
- * "Normalize" a factorial-base number. All of the
- * arithmetic functions call this routine to handle
- * carrys and borrows. A factorial-base number has a
- * proper form where every factorial position in its
- * integer part has a value between 0 and the the
- * magnitude of its position and every factorial position
- * in its fractional part has a value between 0 and one
- * less than the magnitude of its position.
- */
- normalize(n)
- int *n;
- {
- int i, j;
- int *x;
-
- x = &n[HIGHGUARD];
- /*
- * First, check for loss of precision
- * during multiplication or division.
- */
- if(x[LOWGUARD])
- {
- if(nowarning != TRUE)
- fputs("UNDERFLOW\n", stderr);
- x[LOWGUARD] = 0;
- }
- /*
- * Now, work the fractional part first,
- * largest subscript to smallest.
- * The subscript j is the factorial position
- * being put into proper form.
- */
- for(j = MAXFBFRACTION, i = j - 1; i >= 1; --i, --j)
- {
- /* if(x[j] >= j) carry */
- x[i] += (x[j] / j);
- x[j] %= j;
- if(x[j] < 0) /* borrow */
- {
- x[i] -= 1;
- /* modulo= j */
- x[j] += j; /* make positive */
- }
- }
- /* shift any carry to integer part & clear carry */
- n[1] += x[1];
- x[1] = 0;
- /*
- * Now, normalize the integer part,
- * working from smallest subscript to largest.
- * The subscript i is the factorial position
- * being put into proper form.
- */
- x = n;
- for(i = 1, j = 2; i <= MAXFBINTEGER; i++, j++)
- {
- /* if(x[i] >= j) carry */
- x[j] += (x[i] / j);
- x[i] %= j;
- if(x[i] < 0) /* borrow */
- {
- x[j] -= 1;
- x[i] += j;
- }
- }
- if(x[i]) /* if an entry in x[HIGHGUARD] */
- {
- fputs("OVERFLOW\n", stderr);
- x[i] = 0;
- }
- }
-
- /* Add y to x, put result in z */
- add(x, y, z)
- int *x;
- int *y;
- int *z;
- {
- int sign, i;
- int copy[ARRAYSIZE];
-
- if(x[0] != y[0]) /* if different signs */
- {
- if(y[0] == MINUS)
- {
- /*
- * Change the sign of y
- * and subtract y from x.
- */
- negative(y, copy);
- subtract(x, copy, z);
- }
- else
- {
- /*
- * Change the sign of x
- * and subtract x from y.
- */
- negative(x, copy);
- subtract(y, copy, z);
- }
- }
- else
- {
- sign = x[0]; /* save the sign */
- for(i = ARRAYSIZE - 1; i > 0; --i)
- z[i] = x[i] + y[i];
- z[0] = sign;
- normalize(z);
- }
- }
-
- /* Subtract y from x, put result in z */
- subtract(x, y, z)
- int *x;
- int *y;
- int *z;
- {
- int sign, i;
- int copy[ARRAYSIZE];
-
- if(x[0] != y[0]) /* if signs are different */
- {
- negative(y, copy); /* change sign of y */
- sign = x[0]; /* save sign of x */
- add(x, copy, z);
- z[0] = sign;
- return;
- }
- else if(y[0] == MINUS) /* (-x) - (-y) */
- {
- /* if(abs(y) < abs(x)) sign = MINUS */
- sign = lessthan(y, x);
- }
- else
- {
- /* if(x < y) sign = MINUS */
- sign = lessthan(x, y);
- }
- /* Subtract based on the absolute values */
- if(lessthan(x, y))
- {
- for(i = ARRAYSIZE - 1; i > 0; --i)
- z[i] = y[i] - x[i];
- }
- else
- {
- for(i = ARRAYSIZE - 1; i > 0; --i)
- z[i] = x[i] - y[i];
- }
- z[0] = sign;
- normalize(z);
- }
-
- /*
- * Multiply factorial-base x by integer y, result in z.
- * Utility routine called by multiply(), divide(),
- * atofact(), fractoa(), and facttoa(), and always
- * with a positive value for y
- */
- multfbyi(x, y, z)
- int *x;
- int y;
- int *z;
- {
- int i;
-
- for(i = ARRAYSIZE - 1; i > 0; --i)
- z[i] = x[i] * y;
- normalize(z);
- }
-
- /*
- * Divide factorial-base x by integer y, result in z.
- * Utility routine called by multiply(), divide(),
- * and facttoa(), and always with a positive y
- */
- divfbyi(x, y, z)
- int *x;
- int y;
- int *z;
- {
- int i, j, carry, part;
-
- carry = 0;
- /*
- * Work the integer part first,
- * from the largest subscript to smallest
- */
- for(i = MAXFBINTEGER, j = i + 1; i >= 1; --i, --j)
- {
- part = x[i] + carry * j;
- carry = part % y;
- z[i] = part / y;
- }
- /*
- * Now, work the fractional part,
- * from the smallest subscript to largest
- */
- for(i = HIGHGUARD + 1, j = 1;
- j <= MAXFBFRACTION; i++, j++)
- {
- part = x[i] + carry * j;
- carry = part % y;
- z[i] = part / y;
- }
- /*
- * propogate any carry into z[LOWGUARD]
- * to mark underflow and loss of precision
- */
- z[i] = carry * j;
- normalize(z);
- }
-
- /*
- * Multiply factorial-base x by factorial-base y,
- * assign result to z. Uses the identity
- * number * (integer + fraction)
- * == (number * integer) + (number * fraction)
- */
- multiply(x, y, z)
- int *x;
- int *y;
- int *z;
- {
- int i, j, k, sign;
- int partial[ARRAYSIZE];
- int temp[ARRAYSIZE];
- int copy[ARRAYSIZE];
-
- if(x[0] != y[0]) /* if signs different */
- sign = MINUS;
- else
- sign = PLUS;
-
- assignto(partial, zero); /* Initialize result */
- /*
- * Work the integer portion first,
- * from smallest subscript to largest
- */
- absolute(x, copy); /* copy = abs(x) */
- /* first, find largest subscript k where y[k] != 0 */
- for(k = MAXFBINTEGER; (k > 0) && (y[k] == 0); --k)
- ;
- for(i = 1; i <= k; i++)
- {
- /* first shift copy by factorial position */
- multfbyi(copy, i, copy);
- if(y[i]) /* don't bother multiplying by 0 */
- {
- /* multiply by factorial digit */
- multfbyi(copy, y[i], temp);
- add(partial, temp, partial);
- }
- }
-
- /* now work fraction part */
- assignto(copy, x); /* reset copy */
- /* find largest subscript k where y[k] != 0 */
- for(k = ARRAYSIZE - 1; (k > HIGHGUARD + 1) && (y[k] == 0); --k)
- ;
- for(i = HIGHGUARD + 2, j = 2; i <= k; i++, j++)
- {
- /* first shift copy by factorial position */
- divfbyi(copy, j, copy);
- if(y[i]) /* don't bother multiplying by zero */
- {
- /* multiply by factorial digit */
- multfbyi(copy, y[i], temp);
- add(partial, temp, partial);
- }
- }
- partial[0] = sign;
- assignto(z, partial);
- }
-
- /*
- * Divide factorial-base x by factorial-base y, store
- * result in z. Uses blackboard style long division.
- */
- divide(x, y, z)
- int *x;
- int *y;
- int *z;
- {
- int i, j, sign;
- int estimate;
- int copyx[ARRAYSIZE];
- int copyy[ARRAYSIZE];
- int temp[ARRAYSIZE];
- int partial[ARRAYSIZE];
-
- if(x[0] != y[0]) /* if signs different */
- sign = MINUS;
- else
- sign = PLUS;
- absolute(x, copyx);
- absolute(y, copyy);
- assignto(partial, copyx);
- assignto(temp, copyy);
- /*
- * First, estimate the integer part of result by
- * driving y to 1.xxx. Division is VERY slow, so
- * the extra time spent to identify special cases
- * is well worth it.
- */
- if(equalto(temp, zero))
- {
- /* division by zero fault */
- /* not handled here */
- return;
- }
- else if(lessthan(partial, temp))
- {
- assignto(partial, zero); /* integer part 0 */
- }
- else if(lessthan(one, temp))
- {
- /*
- * This could cause a spurious UNDERFLOW
- * message even though the final result
- * would be exact, so we set a flag to
- * suppress the warning.
- */
- nowarning = TRUE;
- while(lessthan(one, temp))
- {
- divfbyi(partial, 2, partial);
- divfbyi(temp, 2, temp);
- }
- multfbyi(partial, 2, partial);
- nowarning = FALSE; /* reset flag */
- }
- else if(lessthan(temp, one))
- {
- while(lessthan(temp, one))
- {
- multfbyi(partial, 2, partial);
- multfbyi(temp, 2, temp);
- }
- }
- else /* division by 1 or -1 */
- {
- assignto(z, x);
- z[0] = sign;
- return;
- }
- /* Now, delete fractional part of estimate */
- for(i = HIGHGUARD + 1; i < ARRAYSIZE; i++)
- partial[i] = 0;
- multiply(copyy, partial, temp);
- while(greaterthan(temp, copyx))
- {
- subtract(partial, one partial);
- subtract(temp, copyy, temp);
- }
- subtract(copyx, temp, copyx);
- multfbyi(copyx, 2, copyx);
- /* partial now holds integer part of result */
-
- /*
- * Now, proceed by long division to divide by
- * the fractional part -- using the subscript
- * (less 1) as the estimate at each position
- */
- for(i = HIGHGUARD + 2, j = 2;
- !equalto(copyx, zero) && (i < ARRAYSIZE);
- i++)
- {
- estimate = j - 1;
- do {
- multfbyi(copyy, estimate--, temp);
- } while(greaterthan(temp, copyx));
- subtract(copyx, temp, copyx);
- partial[i] = ++estimate;
- multfbyi(copyx, ++j, copyx);
- }
- normalize(partial);
- if(sign == MINUS)
- partial[0] = MINUS;
- assignto(z, partial);
- }
-
- /*
- * ASCII to factorial-base conversion
- * Just like ASCII to binary conversion!
- */
- atofact(s, z)
- char s[];
- int *z;
- {
- int i, j, sign;
- int ipart[ARRAYSIZE];
- int fpart[ARRAYSIZE];
-
- assignto(ipart, zero);
- assignto(fpart, zero);
- i = 0;
- sign = PLUS;
- while(s[i] == ' ' || s[i] == '\t')
- i++;
- if(s[i] == '-' || s[i] == '+')
- {
- if(s[i++] == '-')
- sign = MINUS;
- }
- for( ; s[i] >= '0' && s[i] <= '9'; i++)
- {
- multfbyi(ipart, 10, ipart);
- ipart[1] = s[i] - '0';
- normalize(ipart);
- }
- if(s[i] == '.')
- {
- i++;
- j = 0;
- for( ; s[i] >= '0' && s[i] <= '9'; i++)
- {
- multfbyi(fpart, 10, fpart);
- ++j;
- fpart[1] = s[i] - '0';
- normalize(fpart);
- }
- while(j--)
- divfbyi(fpart, 10, fpart);
- add(ipart, fpart, ipart);
- }
- ipart[0] = sign;
- assignto(z, ipart);
- }
-
- /*
- * Convert the fractional part of x to ASCII.
- * Up to count characters go to stdout.
- */
- fractoa(x, count)
- int *x;
- int count;
- {
- int i;
- int temp[ARRAYSIZE];
-
- assignto(temp, x);
- for(i = 1; i <= MAXFBINTEGER; i++)
- temp[i] = 0; /* erase integer part */
- temp[0] = PLUS; /* always positive */
- if(equalto(temp, zero))
- return; /* no fractional part to print out */
- putchar('.');
- while(count-- && !equalto(temp, zero))
- {
- multfbyi(temp, 10, temp);
- putchar('0' + 6 * temp[3] + 2 * temp[2] + temp[1]);
- /* Now erase the integer part */
- temp[3] = temp[2] = temp[1] = 0;
- }
- }
-
- /*
- * CAUTION -- altering the size of outbuff requires
- * some art. If MAXFBINTEGER == 100, it must be
- * large enough to hold the 160 decimal digit integer
- * 9.4259 * 10^159. If MAXFBINTEGER == 180, it must
- * be large enough for the 332 decimal digit integer
- * 3.6362 * 10^331. If you want to deal with really
- * big numbers and increase MAXFBINTEGER, you'll have
- * to give some thought as to how large the conversion
- * buffer is going to have to be.
- */
- /* Allow a little slack */
- char outbuff[MAXFBINTEGER*2] = { 0 };
- int outptr; /* Actually an index and not a "ptr" */
-
- /*
- * Factorial-base to ASCII conversion, integer
- * part and up to count characters of fractional
- * portion go to stdout.
- */
- facttoa(x, count)
- int *x;
- int count;
- {
- int i, j, sign;
- int temp[ARRAYSIZE];
- int val[ARRAYSIZE];
-
- outptr = 0;
- assignto(val, zero);
- assignto(temp, x);
- if((sign = temp[0]) == MINUS)
- {
- temp[0] = PLUS;
- putchar('-');
- }
- for(i = ARRAYSIZE - 1; i > HIGHGUARD; --i)
- temp[i] = 0; /* erase fractional part */
- while(!equalto(temp, zero))
- {
- divbyi(temp, 10, temp);
- for(i = HIGHGUARD + 2, j = 1; j <= 4; i++, j++)
- {
- val[i] = temp[i];
- }
- multfbyi(val, 10, val);
- outbuff[outptr++] =
- '0' + 6 * val[3] + 2 * val[2] + val[1];
- val[3] = val[2] = val[1] = 0;
- /* Now erase fractional portion of temp */
- temp[i-1] = temp[i-2] = temp[i-3] = temp[i-4] = 0;
- }
- if(outptr == 0) /* if no integer part */
- putchar('0');
- else
- {
- while(outptr--)
- putchar(outbuff[outptr]);
- }
- fractoa(x, count); /* to print fractional portion */
- putchar('\n');
- }
-
- /* Remainder of file is RPN calculator & display */
- #ifdef DEBUG
- /*
- * Print a factorial-base number in factorial-base.
- * "digits" are printed between '<' and '>',
- * output goes to stdout.
- */
- facprint(x)
- int *x;
- {
- int i, j;
-
- if(x[0] == MINUS)
- printf("-");
- /* Delete any leading zeroes */
- for(i = MAXFBINTEGER; i >= 1; i--)
- {
- if(x[i] != 0)
- break;
- }
- /* Print any integer portion */
- for( ; i >= 2; )
- printf("<%d>", x[i--]);
- /* Make sure to print at least one digit */
- printf("<%d>", x[1]);
- printf(".");
- /*
- * Print fractional part, deleting any trailing
- * zeroes but printing at least one digit
- */
- i = HIGHGUARD + 2; /* start at 2! */
- printf("<%d>", x[i++]);
- for(j = 0; i < ARRAYSIZE; i++)
- {
- if(x[i] == 0)
- j += 1;
- else
- {
- while(j)
- {
- printf("<0>");
- --j;
- }
- printf("<%d>", x[i]);
- }
- }
- printf("\n");
- }
-
- /*
- * A simple but VERY slow reverse Polish calculator.
- * Commands +, -, *, /, D to print in decimal,
- * F to print in factorial, C to clear the stack,
- * S to display the whole stack in decimal,
- * a decimal number to use as an operand,
- * only 1 operator or operand per line!!
- */
-
- #define IPSIZE 256
- char input[IPSIZE];
- #define STACKSIZE 8
- int stack[STACKSIZE][ARRAYSIZE];
-
- main(argc, argv)
- int argc;
- char *argv[];
- {
- int x, i, prompt, depth;
-
- prompt = (argc > 1) ? TRUE : FALSE;
- depth = 0;
- for( ; ; )
- {
- if(prompt)
- printf(%d> ", depth);
- if(fgets(input, IPSIZE, stdin) == NULL)
- break;
- switch(x = input[0])
- {
- case 'c': /* clear the stack */
- case 'C': depth = 0;
- continue;
- case 'f': /* print top of stack in factorial */
- case 'F': if(depth < 1)
- {
- printf("empty stack\n");
- continue;
- }
- printf("%d: ", depth - 1);
- facprint(&stack[depth-1][0]);
- continue;
- case'd': /* print top of stack in decimal */
- case 'D': if(depth < 1)
- {
- printf("empty stack\n");
- continue;
- }
- printf("%d: ", depth - 1);
- facttoa(&stack[depth-1][0], 50);
- continue;
- case 's': /* display contents of stack */
- case 'S': if(depth < 1)
- {
- printf("empty stack\n");
- continue;
- }
- for(i = 0; i < depth; i++)
- {
- printf("%d: ", i);
- facttoa(&stack[i][0], 50);
- }
- continue;
- case '+': if(depth < 2)
- {
- printf("stack will underflow\n");
- continue;
- }
- add(&stack[depth-2][0],
- &stack[depth-1][0],
- &stack[depth-2][0]);
- --depth;
- continue;
- case '/': if(depth < 2)
- {
- printf("stack will underflow\n");
- continue;
- }
- if(equalto(&stack[depth-1][0], zero))
- {
- printf("division by zero\n");
- /* allow the 0 to be discarded */
- }
- else
- {
- divide(&stack[depth-2][0],
- &stack[depth-1][0],
- &stack[depth-2][0]);
- }
- --depth;
- continue;
- case '*': if(depth < 2)
- {
- printf("stack will underflow\n");
- continue;
- }
- multiply(&stack[depth-2][0],
- &stack[depth-1][0],
- &stack[depth-2][0]);
- --depth;
- continue;
- case '-': if(input[1] == '\n')
- {
- if(depth < 2)
- {
- printf(
- "stack will underflow\n");
- continue;
- }
- subtract(&stack[depth-2][0],
- &stack[depth-1][0],
- &stack[depth-2][0]);
- --depth;
- continue;
- } /* else a negative number */
- else
- break;
- default: if(x != '-' && x != '.' &&
- (x < '0' || x > '9'))
- {
- printf("invalid entry\n");
- continue;
- }
- break;
- /* to convert and stack a number */
- }
- if(depth >= STACKSIZE)
- {
- printf("stack will overflow\n");
- continue;
- }
- atofact(input, &stack[depth++][0]);
- continue;
- }
- }
- #endif
-